/* Updated 10/24/2001. */

/*************************    Program 4.4    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 100

/*void main(){
	
	float matriz[3][3];
	matriz[0][0]= 2;
	matriz[0][1]= 1;
	matriz[0][2]= 1;
	matriz[1][0]= 1;
	matriz[1][1]=1;
	matriz[1][2]=1;
	matriz[2][0]= 2;
	matriz[2][1]=3;
	matriz[2][2]=2;
	
	migs(matriz, 3,);
	
	
}*/


void migs (a,n,x,indx)

/* Function to invert matrix a[][] with the inverse stored
   in x[][] in the output.  Copyright (c) Tao Pang 2001. */
int n;
int indx[];
double a[NMAX][NMAX],x[NMAX][NMAX];
{
  int i,j,k;
  double b[NMAX][NMAX];
  void elgs();

  if (n > NMAX)
  {
    printf("The matrix dimension is too large.\n");
    exit(1);
  }

  for (i = 0; i < n; ++i)
  {
    for (j = 0; j < n; ++j)
    {
      b[i][j] = 0;
    }
  }
  for (i = 0; i < n; ++i)
  {
    b[i][i] = 1;
  }

  elgs (a,n,indx);

  for (i = 0; i < n-1; ++i)
  {
    for (j = i+1; j < n; ++j)
    {
      for (k = 0; k < n; ++k)
      {
        b[indx[j]][k] = b[indx[j]][k]-a[indx[j]][i]*b[indx[i]][k];
      }
    }
  }

for (i = 0; i < n; ++i)
  {
    x[n-1][i] = b[indx[n-1]][i]/a[indx[n-1]][n-1];
    for (j = n-2; j >= 0; j = j-1)
    {
      x[j][i] = b[indx[j]][i];
      for (k = j+1; k < n; ++k)
      {
        x[j][i] = x[j][i]-a[indx[j]][k]*x[k][i];
      }
      x[j][i] = x[j][i]/a[indx[j]][j];
    }
  }
}

void elgs (a,n,indx)

/* Function to perform the partial-pivoting Gaussian elimination.
   a[][] is the original matrix in the input and transformed
   matrix plus the pivoting element ratios below the diagonal
   in the output.  indx[] records the pivoting order.
   Copyright (c) Tao Pang 2001. */

int n;
int indx[];
double a[NMAX][NMAX];
{
  int i, j, k, itmp;
  double c1, pi, pi1, pj;
  double c[NMAX];

  if (n > NMAX)
  {
    printf("The matrix dimension is too large.\n");
    exit(1);
  }

/* Initialize the index */

  for (i = 0; i < n; ++i)
  {
    indx[i] = i;
  }

/* Find the rescaling factors, one from each row */
 
  for (i = 0; i < n; ++i)
  {
    c1 = 0;
    for (j = 0; j < n; ++j)
    {
      if (fabs(a[i][j]) > c1) c1 = fabs(a[i][j]);
    }
    c[i] = c1;
  }

/* Search the pivoting (largest) element from each column */ 

  for (j = 0; j < n-1; ++j)
  {
    pi1 = 0;
    for (i = j; i < n; ++i)
    {
      pi = fabs(a[indx[i]][j])/c[indx[i]];
      if (pi > pi1)
      {
        pi1 = pi;
        k = i;
      }
    }

/* Interchange the rows via indx[] to record pivoting order */

    itmp = indx[j];
    indx[j] = indx[k];
    indx[k] = itmp;
    for (i = j+1; i < n; ++i)
    {
      pj = a[indx[i]][j]/a[indx[j]][j];

/* Record pivoting ratios below the diagonal */

      a[indx[i]][j] = pj;

/* Modify other elements accordingly */

      for (k = j+1; k < n; ++k)
      {
        a[indx[i]][k] = a[indx[i]][k]-pj*a[indx[j]][k];
      }
    }
  }
}
